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An interpretation of AGASA (Akeno Giant Air Shower Array) data by compar- 
ing the experimental results with the simulated ones by CORSIKA (COsmic Ray 
Simulation for KASCADE) has been made. General features of the electromagnetic 
component and low energy muons observed by AGASA can be well reproduced by 
CORSIKA. The form of the lateral distribution of charged particles agrees well with 
the experimental one between a few hundred metres and 2000 m from the core, ir- 
respective of the hadronic interaction model studied and the primary composition 
(proton or iron). It does not depend on the primary energy between 10 17 ' 5 and 
10 20 eV as the experiment shows. If we evaluate the particle density measured by 
scintillators of 5 cm thickness at 600 m from the core (So (600), suffix denotes the 
vertically incident shower) by taking into account the similar conditions as in the 
experiment, the conversion relation from So (600) to the primary energy is expressed 
as -E[eV] = 2.15 x 10 17 So(600) L015 , within 10% uncertainty among the models and 
composition used, which suggests the present AGASA conversion factor is the lower 
limit. Though the form of the muon lateral distribution fits well to the experiment 
within 1000 m from the core, the absolute values change with hadronic interaction 
model and primary composition. The slope of the p^(600) (muon density above 
1 GeV at 600 m from the core) vs. So (600) relation in experiment is flatter than 
that in simulation of any hadronic model and primary composition. Since the ex- 
perimental slope is constant from 10 15 eV to 10 19 eV, we need to study this relation 
in a wide primary energy range to infer the rate of change of chemical composition 
with energy. 
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1 Introduction 



After the observation of cosmic rays with energies greatly exceeding the Greisen- 
Zatsepin-Kuzmin (GZK) cutoff energy [1] by the Fly's Eye [2] and the Akeno Giant 
Air Shower Array (AGASA) [3], seven events exceeding 10 20 eV have been reported 
by the AGASA group [4]. 

To estimate the primary energy of giant air showers observed by the AGASA, the 
particle density at a distance of 600 m from the shower axis (So (600), suffix denotes 
the vertically incident shower) is used, which is known to be a good energy estimator 
[5]. The conversion factor from So(600) [m~ 2 ] to primary energy (E [eV]) at Akeno 
level is derived by simulation [6] based on the COSMOS program by Kasahara et 
at [7] with the QCDjetf 2 ] model [8] and the relation 

E [eV] = 2.0 x 10 17 S (600) L0 (i) 

is used. This relation holds within 20%, independent of primary composition and 
interaction models used [6]. The deviation of the energy spectrum determined in 
this way from that determined by the 1 km 2 array at Akeno is about 10% higher in 
energy around 10 18 eV as shown in [9]. That is, 

£ [eV] = 1.8 x 10 17 S (600) 1 - (2) 

must be applied instead of Eq. 1 to be in agreement with that from the lower energy 
region. 

The AGASA energy spectrum in the highest energy region, which is adjusted to the 
one of the 1 km 2 array at AkencQ [9] are compared with the spectra from the Haverah 
Park [10], the Yakutsk [11] and the stereo Fly's Eye [12] experiments in Fig. 1. All 
four spectra agree with each other within ±15% in energy. It should be noted that 

1 E-mail: mnagano@ccmails.fukui-ut.ac.jp 

2 The QCDjet model is based on QCD and contains a mild scaling violation in the frag- 
mentation region, and a large one in the central region. It includes also production of 
hard mini-jets. This model is extensively used for interpretation of the emulsion chamber 
experiments at Mt.Fuji, Mt. Kambala (China) and Mt. Yanbajin (Tibet) and shows good 
agreement with the experiments in the fragmentation region, where reliable data can not 
be obtained from accelerator experiments. 

3 The spectrum determined between 10 14 ' 5 eV and 10 19 eV using the arrays with differ- 
ent detector spacing at Akeno fits very well with extrapolation of those obtained from 
direct measurement on balloons and satellites, and with the Tibet result obtained through 
the observation of the shower at the height of its maximum development. Therefore the 
AGASA spectrum may be better to normalize to the spectrum by 1 km 2 array. 
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the energy assignment in each experiment has been done separately by each group. 
The Fly's Eye experiment measures the longitudinal development of electrons above 
the observation level and is calorimetric. The Yakutsk group determined the energy 
conversion relation experimentally by measuring not only the lateral distributions of 
electrons and muons but also the energy loss of electrons above the observation level 
from the Cerenkov lateral distribution at the observation level. The Haverah Park 
experiment used water- Cerenkov detectors and the AGASA uses plastic scintillation 
detectors, the energy conversions of both experiments rely on different simulation 
codes. However, the energy assignment is in fairly good agreement as shown by the 
water tank experiment at Akeno [13]. 

On the other hand, Cronin [14] argued that the AGASA energy spectrum should be 
about 30% higher in energy, if the energy is assigned in each experiment according 
to the simulation results based on the MOCCA program with the SIBYLL hadronic 
interaction model. According to their simulation, the relation should be 

E [eV] = 3.0 x 10 17 5 (600) 10 (3) 

at Akeno level. That is, the AGASA energy spectrum results in much higher intensity 
than other experiments at the same energy. 

Recently Kutter pointed out that the difference is due to the large contribution of 
low energy photons [15] in the scintillators in case of AGASA. Though the energy 
loss of electrons and photons in scintillator (p sc ) in units of a vertically traversing 
muon is independent of scintillator thickness, the ratio of p sc to the number of 
charged particles {R sc /ch) depends on the core distance. R sc /ch is about 1.1 within 
100 m from the core and agrees with the experiment [16]. However, it increases as 
the core distance increases. It becomes about 1.2 around 600 m from the core in case 
of MOCCA (SIBYLL) and 1.4 in case of CORSIKA (QGSJET) simulations. The 
difference between the simulation code is due to the difference of energy spectra 
of low energy photons and electrons far from the core. This ratio R sc / C h depends 
also on the number of charged particles which increases as the cutoff energy of 
electromagnetic component decreases. 

Therefore, it is important to evaluate the energy conversion factor Eq. 1 used in 
AGASA experiment with simulation programs other than the COSMOS simulation, 
by taking into account the energy losses of low energy photons and electrons in 
scintillator. 

Since the total number of muons in an extensive air shower with a fixed primary 
energy depends on the chemical composition, the possibility of discriminating heavy 
particles from protons has been extensively studied by many authors. For example, 
Dawson et al. [17] argues that the muon density at 600 m from the core, p^(600), vs. 
energy relation observed at Akeno can be consistently explained by the change of 
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Fig. 1. Energy spectra from the Haverah Park, Yakutsk, stereo Fly's Eye experiments are 
compared with that from the AGASA experiments. 

composition from heavy around 10 17 5 eV to light around 10 19 eV, which is claimed 
by the Fly's Eye group [18]. However, they avoided the discussion of change of 
composition around 10 17 ' 5 eV which is a point of disagreement between the Fly's 
Eye experiment and the Akeno experiment. According to their simulation result 
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[17], the Fly's Eye data shows 100% iron below 10 17 - 5 eV and AGASA data shows 
heavier than iron below 10 17 eV. It is important not only to examine the simulation 
results with other hadronic interaction models, but also to examine the experimental 
relation of p M (600) vs. 5*0(600) in the lower energy region to see whether there is any 
change of slope in the relation. 

As the arrival directions of observed EAS are inclined from the zenith, it is necessary 
to convert 551(600), where 9 represents the zenith angle, to So (600). For zenith angles 
less than 45° (sec# ~1.4), 



50(600) = 50(600) exp 



A 



{sec9 - 1) 



att 



(4) 



is used, where the attenuation length A att = 500 g cm 2 and X is the atmospheric 
depth at Akeno (920 g cnT 2 ) [19]. The attenuation of 5 (600) can be determined 
from integral Se(600) spectra at various 9 by assuming 5*61(600) at constant intensity 
in different zenith angles from primaries of similar energy. If we use data of zenith 
angles up to 55°, 



5 e (600) = 5 (600)exp 



1 n 

— — [sec 9 
Ai v 



1)-|W-1) 2 



(5) 



is better fitted to the data, where Ai = A att = 500 g cm~ 2 and A 2 = 594 g cm~ 2 . 
This was derived, in [20], between 2 x 10 18 eV and 5 x 10 19 eV. Since this attenuation 
with zenith angle is related closely to the properties of hadronic interactions at ultra- 
high energy and to the composition of primary cosmic rays, the interpretation of 
the present result by simulation result is important. 

In this article, we use the EAS simulation program CORSIKA ( COsmic Ray Sim- 
ulation for KASCADE) [21] which was developed at Karlsruhe and is now widely 
distributed and used by various experiments from TeV gamma-rays to the highest 
energy region. Recently the program was improved to be used effectively at the 
highest observed energy [22]. By employing the effective thin sampling procedure, 
the computing time is considerably reduced and hence various combinations of sim- 
ulation conditions with different energies have been realized in this study. 



2 AGASA data 
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2.1 Particle density measured by scintillation detectors on surface 

The AGASA is the Akeno Giant Air Shower Array covering over 100 km 2 area in 
operation at the village of Akeno about 130 km west of Tokyo, to study ultra high 
energy cosmic rays (UHECR) above 10 19 eV [23,24]. 

The lateral distribution of electrons and the shower front structure far from the core, 
which are important to obtain Se(600), are determined with the numerous detectors 
of 1 km 2 array (Al) when showers hit inside the Akeno Branch. The density at core 
distance r is expressed by the function as 

p = N e C e R~ a (l + R)-^-«) (l.O + (^) 2 ) , (6) 

where R = t/Rm, C e is a normalization factor and Rm is a Moliere unit (MU) at 
a height of two radiation lengths above the Akeno level (91.6 m at Akeno) [20]. A 
fixed value of a — 1.2 is used and r\ is expressed as 

77 = 3.97 - 1.79(sec0-l). (7) 

Recently this function was found to be valid up to 3000 m from the core and the 
highest observed energy ~ 10 20 eV [25]. It is important to examine with simula- 
tions whether such an energy independence can be understood with a pure proton 
composition of the primary particles. 

2.2 Muons measured by proportional counters under the absorber 

At AGASA, muons of energies above 0.5 GeV are measured under the lead/iron or 
concrete shielding. At the first stage of AGASA experiment, the lateral distribution 
of muons far from the core was determined with eight muon detectors of 25 m 2 each 
(threshold muon energy : 1 GeV) in the central part of the Akeno Observatory trig- 
gered by the AGASA scintillation detectors on the surface. The lateral distribution 
of muons above 1 GeV is expressed by the following equation [26] . 

ft . = JV,(J)«-™(l + fi)-»(l.0 + (^)*)"'. (8) 

where is a total number of muons and R = r/R Q . is a normalization factor 
and R is a characteristic distance and is expressed by the following equation as a 
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function of zenith angle 6, 

\og(R ) = (0.58 ± 0.04)(sec# - 1) + (2.39 ± 0.05). 



(9) 



This formula can be applied up to 3000 m from the core for showers above 3 x 10 19 eV 
as shown in Doi et al. [27], though the experimental error is still larger than ±50% 
beyond 1000 m from the core in the highest observed energy region. To combine 
the AGASA data of threshold energy of 0.5 GeV with Eq. 8 of 1 GeV threshold, 
the AGASA data is reduced by a factor 1.4, which is determined at 600 m from the 
core, in the whole distance range [28]. 



3 CORSIKA 



In CORSIKA five high energy hadronic interaction models are available and a com- 
parison of available interaction models has been made up to 10 17 eV [29]. Among 
those models QGSJET [30], which includes minijet production in hadronic interac- 
tions, seems to be the best to extrapolate to the highest observed energy range. For 
comparison SIBYLL [31], which is based on a QCD minijet model, is also tried in a 
part of the analysis. 

To see the general aspects of the simulation results in the highest energy region, 
simulations with different interaction models (QGSJET, SIBYLL), primary energies 
(10 17 - 5 , 10 18 - , 10 18 - 5 , 10 19 - , 10 19 - 5 , 10 20 - eV), primary masses (proton, iron), thinning 
levels (10" 5 , 10" 6 ), zenith angles (0°, 29.6° 39.7°, and 51.3° or sec9 = 1.0, 1.15, 1.3, 
and 1.6, respectively), and cutoff energies of electromagnetic component (1.0 MeV, 
0.1 MeV, 0.05 MeV) have been performed. Cutoff energies of hadrons and muons 
are fixed in this series of simulation at 100 MeV and 10 MeV, respectively. 

In each combination, 10 showers are simulated and the average values at the Akeno 
level are summarized. To study fluctuations, 30 showers are used with a thinning 
level of 10~ 6 ' 5 for a limited combinations of conditions. 

In the following all simulation results are from CORSIKA, unless otherwise noted. 



4 Results 
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4-1 Lateral distribution of electrons, photons and muons 

The average lateral distributions of photons, electrons and muons do not depend on 
thinning levels between 10~ 5 and 10~ 6 . However, thinning levels lower than 10~ 6 are 
required to discuss the detailed behavior of scintillator response far from the core, 
taking into account the energy spectrum of individual particles. First the lateral dis- 
tributions of electrons and photons are compared with two different cutoff energies 
of the electromagnetic component, 1.0 MeV and 0.05 MeV. Though there are in- 
creases in the number of electrons and photons with decreasing cutoff energy (about 
1O~2O%0 for the cutoff energies between 1.0 MeV and 0.05 MeV), the change of 
the formQ of the lateral distribution of charged particles is not significant. Therefore 
the cutoff energy of 1.0 MeV is used in the following simulations to save computing 
time and disk space for individual particle information. 



4-1.1 Energy dependence 

In Fig. 2, the primary energy dependence of the lateral distributions of photons, 
electrons and muons is shown for proton primaries with the QGSJET model. Closed 
circles are the lateral distributions of charged particles (muons and electrons) and 
the solid and dashed lines represent the empirical formula of the AG AS A experiment. 
Considering the difference of the primary energy assignment in the experiment and 
the simulation, if we shift the simulated points upward to fit the experimental ones 
within 1000 m from the core, the shapes are in good agreement with experimental 
results up to 1000 ~ 2000 m from the core. 

The energy independence of the lateral distribution of charged particles at the surface 
between 10 17 5 eV and 10 20 eV observed by the AGASA is well supported by the 
QGSJET model. 



4-1.2 Hadronic model dependence 

The lateral distributions of charged particles from SIBYLL are shown in Fig. 3. If 
we shift the simulated values upward, the expected lateral distributions are also 
fitted to the experiment: however, the deviation from experiment slightly increases 
as energy increases (Fig. 3). 



The fraction depends on the primary composition and the hadronic interaction model 
used. 

5 The densities may be shifted vertically to compare the form of the lateral distribution 
with the experimental points, since the primary energy in the experiment is given by the 
conversion equation by Dai et al. 
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Fig. 2. The lateral distributions of photons (open circles), electrons (open squares) and 
muons (open triangles) above 10 MeV and muons (stars) above 1 GeV simulated with 
the QGSJET model. Charged particles which are addition of electrons and muons above 
10 MeV are also plotted by closed circles. The thinning level is 10 -6 and the cutoff energy 
for muons is 10 MeV and those for photons and electrons are 1.0 MeV. Solid and dashed 
lines are the AGASA empirical formula on the surface (Eq. 6) and under absorber (Eq. 
8), respectively. 

While there is no big difference in charged particles for QGSJET and SIBYLL 
calculations, SIBYLL underestimates muons systematically by about a factor 1.6 
with respect to QGSJET, and this would lead to a different interpretation of the 
experimental showers concerning their mass composition. The absolute values are 
related to the assignment of the primary energy and will be discussed in Section 5.3. 
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Fig. 3. Same as Fig. 2 except for the SIBYLL hadronic interaction model being used. 
4-1.3 Primary composition dependence 

In Fig. 4, the lateral distributions of charged particles from iron primary are shown 
for six different primary energies. Energy independence of the form of the lateral 
distribution on primary energy also holds for the iron primary. Therefore the en- 
ergy independence of lateral distribution of AGASA experiment can be understood 
irrespective of primary composition. 

The difference of charged particle densities between proton and iron primaries is 
about 10% around 600 m from the core and those of muons above 1 GeV is about 
50%. 
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Fig. 4. Same as Fig. 2, except for the primary particle being iron. 
4-2 Muon lateral distributions with different cutoff energies 



In Figs. 5 and 6, lateral distributions of muons above 10 MeV, 0.25 GeV, 0.5 GeV, 
1 GeV and 2 GeV are compared with the experimental formula of threshold energy 
of 1 GeV (Eq. 8) at 10 18 eV (left) and 10 19 eV (right). The lateral distributions be- 
come steeper as the cutoff energy increases. While for QGSJET calculations proton 
induced showers describe the experimental curve best, for SIBYLL simulations iron 
induced showers are closer to the experimental distributions. 

It should be noted that the smaller muon number from the SIBYLL model relative 
to other interaction models used in CORSIKA has been pointed out in the 10 14 eV 
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Fig. 5. The lateral distributions of muons of cutoff energies above 10 MeV (top), 0.25 GeV, 
0.5 GeV, 1.0 GeV and 2.0 GeV (bottom). QGSJET model. The thinning level is 10~ 6 . 
Upper figures are proton primaries and lower ones are iron. The incident energies are 10 18 
eV (left) and 10 19 eV (right), respectively. Symbols as in Fig. 2. 

~ 10 15 eV region by Knapp et al. [29]. The slope of energy spectra of muons between 
0.25 GeV and 1.5 GeV at 600 m from the core is similar irrespective of the QGSJET 
and the SIBYLL models and proton and iron primaries for 10 19 eV. The experimental 
result [28] agrees well with the simulated ones except at 0.25 GeV, where we have 
in the experiment punch-through of the electromagnetic component into the muon 
detectors, even at 600 m from the core. 

The relation of muon density above 1 GeV at 600 m from the core (p At (600)) with the 
primary energy for three different combinations of simulation codes and interaction 
models is compared in Table 1. Here the slope a in p M oc E a and the density at 
10 19 eV are listed. The results from the SIBYLL model with the MOCCA code [17] 
are also listed for comparison. Though the hadronic model is the same, there is about 
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Fig. 6. Same as Fig. 5, except for the hadronic interaction model is being SIBYLL. 

10% difference in density between the results based on the two different simulation 
codes. The implications from this comparison will be made in Section 5.4. 



4-3 Charged particle density, Sq(600) 



The relation of 5*0(600) with the primary energy is examined for various conditions. 
5*0(600) does not depend on interaction model (QGSJET, SIBYLL) nor primary 
mass (proton, iron) within 20% which supports the previous simulations by Hillas 
et al. [5] and Dai et al. [6]. The 5o(600) vs. primary energy relation is almost linear, 
irrespective of primary mass (proton or iron) in case of QGSJET. However, 5o(600) 
~ E ' 97 in case of SIBYLL. There are some differences between the simulated number 
of charged particles and the relation used at AGASA. The details will be discussed 
in Section 5.3. 
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Table 1 

Comparison of slope and density at 10 19 eV 



Code 


Model 


Primary 


Slope 


p M (600) 
at 10 19 eV 


Note 


CORSIKA 


QGSJET 


proton 


0.92 


3.85 




CORSIKA 


QGSJET 


iron 


0.89 


5.64 




CORSIKA 


SIBYLL 


proton 


0.88 


2.39 




CORSIKA 


SIBYLL 


iron 


0.87 


3.96 




MOCCA 


SIBYLL 


proton 


0.90 


2.95 


[17] 


MOCCA 


SIBYLL 


iron 


0.86 


4.57 


[17] 



4-4 Zenith angle dependence of charged particles 

There is a zenith angle dependence of the lateral distribution of charged particles 
in the AGASA experiment and the slope parameter r\ is expressed by Eq. 7 as a 
function of zenith angle. In Fig. 7, the lateral distribution of charged particles at a 
zenith angle of 39.7° is shown for proton primaries of six primary energies. For the 
AGASA lateral distribution, rj = 3.6 and an attenuation length A att = 500 g cm -2 
are used. It is shown that the simulated results fit well to the experimental function, 
including the absolute values, independent of primary energy. 

Figure 8 shows the lateral distributions of charged particles at the four zenith angles 
0°, 29.6°, 39.7° and 51.3° for proton (left) and for iron (right) primaries and for 
primary energies 10 18 ' 5 eV and 10 19 ' 5 eV, respectively. It should be noted that there 
are no differences in attenuation of charged particles with zenith angle between 
proton and iron primaries. The details will be discussed in Section 5.5. 



5 Discussion 

The general features of the electromagnetic component and the low energy muons 
observed by AGASA can be well reproduced by CORSIKA simulation, except the 
slight differences in the absolute values. These discrepancies in absolute values be- 
tween experiment and simulation are partly due to the assignment of primary energy, 
which depends on the definition of a single particle used in experiment and simu- 
lation. More remarkable discrepancies are the slope of p M vs. 5*0(600) relation and 
zenith angle dependence 5*0(600) for constant primary energy, which are shown in 
Fig. 12 and in Fig. 14, respectively. In the following we evaluate the simulated results 
with a similar definition of each observable as in the AGASA experiment. 
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Fig. 7. The lateral distributions of photons, electrons, muons and charged particles simu- 
lated at a zenith angle of 39.7°. The thinning level is 1CT 6 and the cutoff energy of muons 
is 10 MeV and those of electrons and photons are 1 MeV. For drawing the AGASA lateral 
distribution, an attenuation length of 500 g cm -2 is used and the change of n with zenith 
angle is taken into account. Symbols as in Fig. 2. 



5.1 Definition of density used in Akeno/AGASA experiment 



In the AGASA experiment, a scintillator of 5 cm thickness is used to detect particles 
on the surface. The scintillator is placed inside an enclosure made of iron of 1.5 ~ 
2 mm thickness and the detector is in a hut whose roof is made of an iron plate of 
0.4 mm thickness. 
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Fig. 8. The lateral distributions of charged particles at four zenith angles 0°, 29.6°, 39.7°, 
and 51.3° for primary energies 10 18 ' 5 eV (top) and 10 19 ' 5 eV (bottom). QGSJET model 
and cutoff energy of muons is 10 MeV and those of photons and electrons are 1.0 MeV. 
The left figure is for proton primary and the right one is for iron primary. 



The definition of a single particle at the Akeno experiment (Vi) is the average of 
the pulse- height distribution of muons traversing a scintillator vertically [9]. Vi is 
10% larger than the peak value, since the distribution is not Gaussian, but subject 
to Landau fluctuations. If we use the peak value in pulse height distribution (PHD) 
of omnidirectional particles, the peak value V p h is accidentally coincident with Vi 
at Akeno level (V p h — Vi). To convert a density measured by a scintillator to an 
electron density corresponding to the calculated density by the NKG function [32], 
the density measured by scintillators and spark chambers at the Institute for Nuclear 
Study at Tokyo [16] is used at Akeno. This ratio is 1.1 between 10 m and 100 m 
from the core and hence the electron shower size (N e ) determined by the Akeno 1 
km 2 array was reduced by 10% from the calculated N e [9]. 
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In AGASA a peak value V pw in the pulse width distribution (PWD) of omnidirec- 
tional muons on a 5 cm scintillator is used as a single particle conventionally. The 
pulse width is obtained by discriminating a signal, which decays exponentially with 
a time constant of 10 //sec, at a constant level. V pw is not equal to V p h and is related 
to V p h as : 

V pw = (V ph + ^V p l + ^)/2 

where a is a full width at half maximum in PHD [33]. By putting V p h = 1.0 and a = 
0.70, we obtain V pw = 1.1. A conventional value, V pw , used at AGASA corresponds 
to l.lxVi(~ 1.1 x V^). That is V pw corresponds to a measured density by spark 
chamber or the electron density, as far as the ratio of the density measured by 
scintillators and the spark chambers is 1.1. 

Though there is a transition effect of the electromagnetic component in scintillator 
of 5 cm thickness within 30 m from the core [9], the density of charged particles as 
expressed in units of Vi does not depend on the thickness of the scintillator above 
30 m from the core as shown experimentally in Teshima et al. [34]. This can be 
understood since the radiation lengths of scintillator and air are very similar, so 
that the fraction of electrons at any depth in the scintillator changes only slowly as 
compared to air. This independence of thickness of scintillator has also been shown 
in the simulations of Cronin [14] and Kutter [15]. 

Assuming the 10% difference of scintillator density to spark chamber density at 600 
m from the core, the density in units of V pw (1.1 x V p h) coincide with spark chamber 
density as described above and Eq. 1 was applied to deriving the AGASA energy 
spectrum. 

5.2 Densities measured by scintillator of 5cm thickness 

So far the AGASA group has used a factor 1.1 of scintillator density to spark chamber 
density which is determined within 100m from the core [16]. However, the factor has 
not yet been measured beyond 100 m from the core. In the following we discuss 
the lateral distribution of energy losses by photons, electrons and muons in 5cm 
scintillator in units of V pw at Akeno level, taking into account the incident angles 
of electrons and photons far from the core, because the incident angles of electrons 
and photons on scintillator may not be vertical even for a vertical shower and the 
particles have some angular distribution. 

In Fig. 9 the energy spectra of photons, electrons and muons are shown at core 
distances between 500 m and 800 m, simulated by CORSIKA. Many photons still 
remain at a zenith angle of 51.3° (sec# = 1.6). 
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Fig. 9. Energy spectra of photons, electrons and muons at core distance between 500 m 
and 800 m. The primary particle is a proton of 10 19 eV. Zenith angles are 0° (left) and 
51° (right). 



Using CORSIKA, in a shower the incident angle of each particle to the surface is 
recorded, so that we can calculate the energy loss of each particle which is incident 
on the scintillator with various zenith angles. In the case of electrons and muons, 
only ionization loss is taken into account. The energy loss of photons is evaluated 
as follows. By dividing the scintillator in thin layers, the fraction of conversion to 
electrons in each layer is calculated using the photon attenuation length in water 
[35]. The energy loss of electron or electron/positron in the remaining layers for 
the fraction of photons is calculated. In this way the average energy loss of a single 
photon in a scintillator of various thicknesses is calculated as a function of the photon 
energy. Though the calculation process is simple, the result for a scintillator of 5 cm 
thickness agrees well with the Monte Carlo simulation results by Kutter [15]. The 
density of a shower of zenith angle 9 is evaluated as the energy loss in scintillator of 
5 x sec# cm thickness (in units of V pw ) in an area of 1 x cos# m 2 . 

In Fig. 10, the lateral distribution of energy deposit in scintillator of 5 cm thickness 
in units of V pw (p sc ) is plotted by closed circles and compared with that of the 
experimental lateral distribution (dashed line). The simulated lateral distribution is 
flatter than the experimental one. p sc reflects the number of electrons near the core 
(up to about 200 m from the core), but becomes larger than the electron density 
with core distance. Though the agreement of p sc (600) with the experimental So(600) 
is quite good, the difference of the lateral distribution of p sc from the experiment 
must be studied further. 
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Fig. 10. Lateral distribution of energy deposit in scintillator (•) of 5 cm thickness (p sc ) in 
units of V pw is compared with the experimental lateral distribution of AGASA. Those of 
electrons (> 1 MeV), photons (> 1 MeV), muons (> 0.5 GeV) and muons (> 1 GeV) are 
also shown. The left figure is proton primary and the right one iron. 

5.3 Primary energy and Sq(600) relation evaluated by CORSIKA 

In Table 2, the density of charged particles at 600 m from the core, p c / l (600), or the 
scintillator density of 5 cm thickness, p sc (600), are listed for showers of 10 19 eV of 
vertical incidence. 

The various combinations of simulation codes (CORSIKA, MOCCA), hadronic inter- 
action models (QGSJET, SIBYLL), primary species (proton, iron), thinning levels 
and threshold energies of the electromagnetic components are compared. 

In general, the difference of p c / l (600) due to the difference of simulation codes or 
hadronic interaction models is within 10% for the same cutoff energy of electromag- 
netic component. p c h(600) depends on the cutoff energy of electromagnetic compo- 
nent. In the 5*o (600) vs. energy relation by Dai et al. (Eq. 1), the electromagnetic 
component with energy of less than the thinning level is connected to the NKG 
function without cutoff energy for electrons and photons, and MU at 2 radiation 
lengths above the observation level is used. The result by CORSIKA simulated with 
the similar method is given in Table 2 as QGSJET-NKG. In case of CORSIKA, the 
relation is 

E[eV] = 2.2 x lO^So^OO) 1 - (10) 
and is about 10% larger than that by Dai et al. 

The energy losses of photons, electrons and muons in scintillator of 5 cm thickness 
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Table 2 

Comparison of charged particle density (p c h) or scintillator density in units of Vi (p sc ) at 
600 m from the core for showers of vertical incidence, 10 19 eV and proton or iron primary. 



Code 


"TV J IT 

Model 


Primary 


mi ■ 

1 running 


mi i 11 

I hreshold 


p c/l (600) or 


Note 








1 i 

level 


E e7 (MeV) 


p sc (600) (m z ) 




CORSIKA 


QGSJET 


proton 


10~ 5 


1.0 


32.5 




CORSIKA 


QGSJET 


proton 


10~ 6 


1.0 


32.4 




CORSIKA 


QGSJET 


proton 


10~ 6 


0.1 


37.5 




CORSIKA 


QGSJET 


proton 


10~ 6 


0.05 


39.1 




CORSIKA 


QGSJET 


iron 


10~ 5 


1.0 


35.7 




CORSIKA 


QGSJET 


iron 


10~ 5 


0.05 


39.8 




CORSIKA 


SIBYLL 


proton 


10~ 6 


1.0 


30.4 




CORSIKA 


SIBYLL 


iron 


10~ 6 


1.0 


33.5 




CORSIKA 


QGSJET 


proton 


10~ 6 


1.0 


38.0 















45.0 




UUKolrvA 


Ann TTTT^ 

QLoJrijl 


proton 


10 


1.0 


«i /ion 

sci. 43.0 




V V 7 I V . ' J 1 \ . \ 




iron 


i n -5 


1 n 


sci. ^u.z 






QTT5VT T 


proton 


i c\— 6 
1U 


1 n 
l.U 


SCI. OO.Z 








iron 




1 n 

l.U 


Q „; AAA 
SCI. 44.4 




MOCCA 


SIBYLL 


proton 




0.1 


33.5 


Cronin W 






iron 




0.1 


38.7 


Cronin ^ 


COSMOS 


QCDjet 


proton 


10~ 5 





50.0 


Dai et al.( 2 ) 




-NKG 


CNO 


10~ 5 













iron 


10~ 5 










W Simulation results made at Fermi Lab. using the SIBYLL interaction model with 
MOCCA simulation code (J. Cronin [14]). 



( 2 ) Two dimensional simulation results made at ICRR with COSMOS by Dai et al. [6]. 
Photons and electrons of energies below the thinning energy level are connected to the 
NKG function in which the Moliere length is used at 2 radiation lengths above the Akeno 
altitude. 

(p sc (600)) in units of V pw have been evaluated as described in the previous section 
by taking into account their incident angle to the scintillator and attenuation of low 
energy photons and electrons in the scintillator container and hut. The relation is 
drawn in Fig. 11 for proton and iron primaries. For the QGSJET hadronic interaction 
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D Proton QGSJET 
A Iron QGSJET 

O Proton SlBYLL 
Iron SlBYLL 



19.5 20 20.5 

log(ENERGY eV) 

Fig. 11. 5o(600) or p sc (600) vs. energy relation estimated from the QGSJET and SlBYLL 
model for proton and iron primary. 



model we obtain 



E[eV] = 2.07 x 10 17 p sc (600) 103 for proton (11) 

E[eV] =2.24 x 10 17 p sc (600) 100 for iron (12) 

and for SlBYLL the relation reads 

E[eV] =2.30 x 10 17 p sc (600) 103 for proton (13) 

E[eV] =2.19 x 10 17 p sc (600) 101 for iron. (14) 



Though there is a difference between proton and iron showers or QGSJET and 
SlBYLL, any combination assigns a higher energy than that by Eq. 1. The functions 
given by Eqs. 11, 12 and 14 intersect around p sc (600) = 20 ~ 30, the combined 
relation may be written as 



E[eV] = 5.65xl0 18 f^#^V' 015 , (15) 



25.0 

by taking the average slope of the four equations. If we rewrite the above equation 
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in similar form as Eq. 1, the following function may be used. 

E[eV] = 2.15 x 10 17 S (600) L015 . (16) 

Here we denoted So (600) instead of p sc (600) as a scintillator density used in the 
experiment, since S (600) ~ p sc (600) as described before. This equation means that 
the energy of the AGASA events must be increased by 14% at 10 19 eV and 18% at 
10 20 eV, if we evaluate the primary energy by CORSIKA. 

5.4 Composition deduced from the muon component 

The slope of the p^(600) vs. energy relation measured by experiment is smaller 
than that in CORSIKA simulations. In Fig. 12 the p M (600) vs. S (600) (or p sc (600)) 
relation is shown. Since the relation So(600) vs. energy is almost linear irrespective 
of primary composition or interaction model used, the relation 

p^(600) = A x S (600) b 

is similar to p^(600) vs. energy relation and the parameter b is similar to a in Table 
1. 

To explain the difference of slopes between experiment and simulation, Dawson et 
al. [17] favours a likely change in mass composition from heavy to light at energies 
above 10 18 eV, as a supporting evidence of the change of mass composition claimed 
by the Fly's Eye experiment [18]. The difficulty of their interpretation arises when 
the relation is applied to the energy region lower than 10 17 eVQ. The iV M vs N e 
relation of the Akeno result [26] is compared with the result of KASCADE [36] in 
Fig. 13 without any correction for each experiment. The empirical formula of Akeno 
at sec^ = 1.05 is expressed by 

= (2.94 ± 0.14) x 10 5 x (iV e /10 7 ) - 76±0 - 02 . (17) 

N e attenuates about a factor 1.5 from the Akeno level (920 g cm -2 ) to the KASCADE 
level (1020 g cm" 2 ), while of the KASCADE experiment with the threshold 

6 The AGASA group claims that no change of composition around 10 17 ' 5 eV is observed 
by their experiment [26] in contrast to the Fly's Eye experiment. Dawson et al. do not 
argue on this point (compare with Fig. 4 and 6 of their paper). According to their figure, 
the fraction of iron component is 100% below 10 17 ' 5 eV in the Fly's Eye experiment, while 
that is 80% at 10 17 5 eV and exceeds 100% below 10 17 - eV in the AGASA experiment. 
There is still no agreement between the two experiments below 10 17 ' 5 eV, which is the 
point in dispute. 
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Fig. 12. Comparison of experimental vs So (600) or p c /j relation with that from COR- 
SIKA simulation. 

energies of 2 GeV is expected to be smaller than that at Akeno of 1 GeV about 
a factor 1.2 ~ 1.5 [37]. Therefore such an agreement of the extrapolation of the 
absolute values from both experiments is well understood. The important result, 
however, is the agreement of the slopes of the vs. N e relation of both experiments 
in quite different energy regions. 

Since in the higher energy region, N e can not be determined by AGASA, the p M (600) 
vs. 5*0(600) relation is used. The result at sec^ = 1.09 [26] is expressed as 

p^(600) = (0.16 ± 0.01) x 5 (600) a82±a03 . (18) 

Taking into account the relation S , (600) oc E oc A^?- 9 , the above equation coincides 
with Eq. 17 in the overlapping energy region [26]. Therefore the slope seems not to 
change from 10 14 ' 5 eV to 10 19 eV. If we take the absolute values of SIBYLL model as 
used in Dawson et al. [17] (refer to Fig. 4 of their paper), the composition becomes 
heavier than iron below 10 17 eV. This conclusion demonstrates that it is important 
to compare the experimental results with simulations in as wide energy range as 
possible. 



24 



^ 8 

* 7 
6.5 

6 
5.5 

5 
4.5 

4 
3.5 



-•- Akeno empirical relation 

• Akeno 1 km 2 array > 1 GeV 

□ KASCADE E^ > 400 MeV 

■ KASCADE E . > 2 GeV 



□ 

° ■ 



□ j 
□ ■' 



□ i 
□ »- 



□ ■ 
□ ■ 



4 5 6 7 8 9 10 

log(N e ) 

Fig. 13. vs. N e relation of the Akeno and the KASCADE experiments. The solid line 
represents Eq. 16 and the dashed one is its extrapolation. 

5.5 Attenuation of Sq(600) with zenith angle and the implication on the primary 
composition around 10 19 eV 



The p sc (600) values are plotted in Fig. 14 as a function of atmospheric depth (sec#) 
and connected with lines. The solid lines represent proton primaries and dotted ones 
are iron primaries. The experimental points of AGASA [38] are also plotted. These 
were obtained using the method of 'equi-intensity cuts' on the integral p sc (600) 
spectra, based on the assumption that the flux of showers above a certain primary 
energy does not change with atmospheric depth. As shown in the figure the variation 
of p sc (600) with zenith angle by CORSIKA (QGSJET) is similar to proton and iron 
primaries and is smaller than that of the experiment. If we take into account the 
error in p ch (Q00) and zenith angle determination, this difference becomes larger, since 
flux of the observed p c h (600) spectrum is increased by a factor of exp(<r 2 (7 — l) 2 /2), 
where a is the error in p c h(600) determination in a logarithmic scale and 7 is the 
power index of the differential p c /i (600) spectrum [39]. 

To examine the difference, the zenith angle variation of the electron density (> 1 
MeV) at 600 m from the core and the muon density (> 1 GeV) at 600 m from 
the core are plotted in Fig. 15. The irregularities in the curves may be due to the 
statistical fluctuations of the limited number of simulated showers at each point 
and fitting errors to derive the density at 600 m. Within these uncertainties the 
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Fig. 14. The variation of p sc (600) with zenith angle. Solid lines are from protons and dotted 
are from iron (QGSJET model, thinning level 10~ 6 ). Experimental points [20] correspond 
to 40, 20, 10, 5 and 2 x 10 18 eV from top to bottom, respectively. 

attenuation with zenith angle by CORSIKA and that by Dai et al. agree with each 
other as in the left figure of Fig. 15. The difference of absolute values is partly due 
to the difference of cutoff energies as described before. Since the number of muons 
does not change with zenith angle and the number of muons in iron initiated showers 
are larger than that in proton showers as shown in the right figure of Fig. 15, the 
attenuation of p c / l (600), which consists of electrons above 1 MeV and muons above 
10 MeV, does not depend on primary energy and composition. 

As is described before, 5*0(600) must be evaluated as p sc (600), energy deposit in a 
scintillator in units of V pw , at various zenith angles taking into account the increase 
of scintillator thickness with zenith angle. It is found that the attenuation of elec- 
tromagnetic density with zenith angle is compensated with muon density and hence 
the attenuation of p sc (600) with zenith angle is similar for proton and iron primaries 
in case of the QGSJET model. Therefore the muon density relative to the electro- 
magnetic density, which is sensitive to the hadronic model, is an important quantity 
to compare with the experiment. Since the AGASA data have been increased con- 
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Fig. 15. The variation of electron density at 600 m from the core with zenith angle (left) 
and that of muon density above 1 GeV at 600 m from the core (right). Solid lines are 
from proton primaries and dotted are from iron (QGSJET model, thinning level 10 -6 and 
cutoff energy of electrons are 1.0 MeV). Triangles connected with dotted lines are from 
Dai et al. [6] for primary protons of 10 18 ' eV and 10 190 eV, respectively. 

siderably at the highest energy region compared to the one used in the Fig. 14, it 
may be better to await further comparisons for the publication of new experimental 
results. 



6 Conclusion 



An interpretation of AGASA data has been made by comparing them with results 
simulated with CORSIKA. Ten air showers for each of the primary energies 10 17 5 
eV, 10 18 -° eV 10 18 - 5 eV, 10 19 eV 10 19 ' 5 eV and 10 20 eV are simulated under each 
combination of the QGSJET or SIBYLL hadronic models and proton or iron pri- 
maries. General features of the electromagnetic component and low energy muons 
observed by AGASA can be well reproduced by CORSIKA. There are some discrep- 
ancies which must be studied in more detail to improve the agreement between the 
experiment and simulation. In the following some results implicated by the simula- 
tion related to the AGASA experiment are summarized. 

The form of the lateral distribution of charged particles agrees well with the exper- 
imental one between a few hundred m and 2000 m from the core, irrespective of 
hadronic interaction model or composition and it does not depend on the primary 
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energy simulated. Though the shape of the muon lateral distribution fits also the ex- 
periment, absolute values change with the hadronic interaction model and primary 
composition. 

If we evaluate the density measured by scintillator of 5 cm thickness as 5*0(600) by 
taking into account similar conditions as in the experiment, the conversion relation 
from 50(600) to the primary energy can be written as 

E[eV] = 2.15 x 10 17 S (600) L015 . 

This means that the energy assignment of the AGASA experiment is shifted to 
higher energies, if we estimate it using CORSIKA. 

The slope of p M vs. So (600) relation in the experiment is flatter than that in simu- 
lations of any hadronic model and primary composition. The situation may not be 
changed even by taking into account the primary energy spectrum and fluctuations 
of p M and So (600). Since the slope seems to be constant in a wide primary energy 
range, we need to study this relation over a wide energy range. Otherwise the com- 
position may be interpreted as heavier than iron or lighter than proton outside the 
narrow investigated energy region. 

There is a disagreement of the attenuation length determined at AGASA and that 
by CORSIKA simulation, even if we take into account the particle incident angle 
to the scintillator. If we take into account the experimental error in zenith angle 
determination and So (600) determination, the disagreement increases. Since the ex- 
perimental data used in the present analysis is still preliminary, we better wait for 
the publication of new experimental values for further discussions. 
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